TP4 - Clustering
library(FactoMineR)
library(ggplot2)
library(gridExtra)
library(factoextra)
library(cluster)
library(mclust)
library(fpc)
library(phyloseq)Description des données
Les données étudiées dans ce TP sont issues d’analyses chimiques de plusieurs vins de la même région de l’Italie mais correspondant à trois cépages différents. Les analyses chimiques ont permis de mesurer les quantités de 13 constituants contenus dans chacun de ces 3 types de vins.
Question 1
Commencez par lire les données et vérifier si la lecture s’est bien passée à l’aide des commandes suivantes.
url <- url("https://raw.githubusercontent.com/cmaugis/Stat-Avancees---PFB-Toulouse/main/Data/WineTP.txt")
Wine = read.table(url)
Wine[, 1] <- as.factor(Wine$label)
dim(Wine)[1] 178 14
label Alcohol Malic.acid Ash Alcalinity.of.ash Magnesium Total.phenols
1 1 14.23 1.71 2.43 15.6 127 2.80
2 1 13.20 1.78 2.14 11.2 100 2.65
3 1 13.16 2.36 2.67 18.6 101 2.80
4 1 14.37 1.95 2.50 16.8 113 3.85
5 1 13.24 2.59 2.87 21.0 118 2.80
Flavanoids Nonflavanoid.phenols Proanthocyanins Color.intensity Hue
1 3.06 0.28 2.29 5.64 1.04
2 2.76 0.26 1.28 4.38 1.05
3 3.24 0.30 2.81 5.68 1.03
4 3.49 0.24 2.18 7.80 0.86
5 2.69 0.39 1.82 4.32 1.04
OD280.OD315 Proline
1 3.92 1065
2 3.40 1050
3 3.17 1185
4 3.45 1480
5 2.93 735
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
Question 2
Enoncé
A l’aide de la fonction table() déterminez la répartition des 178 vins en 3 cépages.
Question 3
Enoncé
A l’aide des commandes suivantes, faites une ACP des données. Commentez.
Correction
ACP des données avec couleur selon les labels :
respca <- PCA(Wine, scale.unit = TRUE, quali.sup = 1, graph = FALSE)
plot(respca, choix = "ind", habillage = 1, label = "none")Commentaire : On repère bien les 3 cépages dans le premier plan factoriel. Grâce au deuxième graphique, on peut regarder comment les variables intiales sont corrélées aux deux premières composantes principales. On oublie dans la suite cette information pour faire la classification non supervisée des données et on se comparera avec la variable cépage.
Méthodes par partitionnement
Question 1
Enoncé
A l’aide de la fonction kmeans(), déterminez une classification en \(K=3\) groupes avec la méthode des Kmeans.
Vous pouvez visualiser la classification obtenue avec la commande suivante :
A l’aide des commandes suivantes, comparez la classification obtenue avec la répartition en cépages.
Correction
Classification en \(K=3\) groupes avec la méthode des Kmeans :
Length Class Mode
cluster 178 -none- numeric
centers 39 -none- numeric
totss 1 -none- numeric
withinss 3 -none- numeric
tot.withinss 1 -none- numeric
betweenss 1 -none- numeric
size 3 -none- numeric
iter 1 -none- numeric
ifault 1 -none- numeric
Visualisation de la classification obtenue :
1 2 3
1 13 0 46
2 20 50 1
3 29 19 0
[1] 0.3711137
[1] 0.2977528
Question 2
Enoncé
A l’aide de la fonction pam(), déterminez une classification en \(K=3\) groupes avec la méthode PAM.
Visualisez la classification obtenue et comparez-la avec la répartition en cépage.
Comparez cette classification avec celle des Kmeans :
Correction
On cherche une classification en \(K=3\) classes avec PAM :
Les médoides des \(3\) classes sont :
Alcohol Malic.acid Ash Alcalinity.of.ash Magnesium Total.phenols
51 13.05 1.73 2.04 12.4 92 2.72
136 12.60 2.46 2.20 18.5 94 1.62
73 13.49 1.66 2.24 24.0 87 1.88
Flavanoids Nonflavanoid.phenols Proanthocyanins Color.intensity Hue
51 3.27 0.17 2.91 7.20 1.12
136 0.66 0.63 0.94 7.10 0.73
73 1.84 0.27 1.03 3.74 0.98
OD280.OD315 Proline
51 2.91 1150
136 1.58 695
73 2.78 472
[1] 51 136 73
On peut visualiser cette classification dans le premier plan factoriel :
On compare cette classification avec les cépages et avec celle des Kmeans :
[1] 0.3715001
1 2 3
1 46 2 0
2 13 19 30
3 0 50 18
[1] 0.9663286
1 2 3
1 1 0 47
2 61 1 0
3 0 68 0
Question 3
Enoncé
Dans cette question, on souhaite sélectionner le nombre de classes pour la méthode des Kmeans. A l’aide des commandes suivantes, comparez les différents critères.
PlotCriteres <- function(Data, Kmax) {
Iintra = NULL
CH = NULL
Silhou = NULL
for (k in 2:Kmax) {
Classif = kmeans(Data, k, nstart = 10)
Iintra = c(Iintra, sum(Classif$withinss)/nrow(Data))
aux = silhouette(Classif$cl, daisy(Data))
Silhou = c(Silhou, mean(aux[, 3]))
CH = c(CH, calinhara(Data, Classif$cluster))
}
dfindice = data.frame(K = 2:Kmax, CH = CH, Iintra = Iintra, Silhouette = Silhou)
g1 = ggplot(dfindice, aes(x = K, y = Iintra)) + geom_point() + geom_line() +
ggtitle("Iintra") + theme(plot.title = element_text(size = 9)) + ylab("")
g2 = ggplot(dfindice, aes(x = K, y = CH)) + geom_point() + geom_line() + ggtitle("Calinski-Harabasz") +
theme(plot.title = element_text(size = 9)) + ylab("")
g3 = ggplot(dfindice, aes(x = K, y = Silhouette)) + geom_point() + geom_line() +
ggtitle("Silhouette") + theme(plot.title = element_text(size = 9)) + ylab("")
grid.arrange(g1, g2, g3, ncol = 3)
}Correction
PlotCriteres <- function(Data, Kmax) {
Iintra = NULL
CH = NULL
Silhou = NULL
for (k in 2:Kmax) {
Classif = kmeans(Data, k, nstart = 10)
Iintra = c(Iintra, sum(Classif$withinss)/nrow(Data))
aux = silhouette(Classif$cl, daisy(Data))
Silhou = c(Silhou, mean(aux[, 3]))
CH = c(CH, calinhara(Data, Classif$cluster))
}
dfindice = data.frame(K = 2:Kmax, CH = CH, Iintra = Iintra, Silhouette = Silhou)
g1 = ggplot(dfindice, aes(x = K, y = Iintra)) + geom_point() + geom_line() +
ggtitle("Iintra") + theme(plot.title = element_text(size = 9)) + ylab("")
g2 = ggplot(dfindice, aes(x = K, y = CH)) + geom_point() + geom_line() + ggtitle("Calinski-Harabasz") +
theme(plot.title = element_text(size = 9)) + ylab("")
g3 = ggplot(dfindice, aes(x = K, y = Silhouette)) + geom_point() + geom_line() +
ggtitle("Silhouette") + theme(plot.title = element_text(size = 9)) + ylab("")
grid.arrange(g1, g2, g3, ncol = 3)
}gskmn <- clusGap(Wine[, -1], FUN = kmeans, nstart = 10, K.max = 10, B = 60)
plot_clusgap(gskmn) + theme(plot.title = element_text(size = 9)) + xlab("K")Avec l’inertie intra, on peut choisir 4 ou 7 classes. Avec Silhouette, on choisit 2 classes.
Question 4
Enoncé
Etudiez la classification retenue avec le critère Silhouette.
Question 5
Enoncé
En adaptant des commandes des questions précédentes, déterminez une classification des vins à partir des coordonnées de l’ACP à l’aide de la méthode des Kmeans.
Correction
gskmn <- clusGap(DataACP, FUN = kmeans, nstart = 10, K.max = 10, B = 60)
plot_clusgap(gskmn) + theme(plot.title = element_text(size = 9)) + xlab("K")ClassifKmeansPCA = kmeans(DataACP, 3)
fviz_cluster(ClassifKmeansPCA, data = Wine[, -1], ellipse.type = "norm", labelsize = 8) +
ggtitle("")
1 2 3
1 59 0 0
2 3 65 3
3 0 0 48
[1] 0.897495
Classification hiérarchique
Question 1
Enoncé
Comparez les classifications en 3 classes obtenues avec une CAH avec le lien de ward etle lien complet respectivement.
Correction
Kfixe = 3
d <- dist(Wine[, -1], method = "euclidian")
hward <- hclust(d, method = "ward.D2")
clward <- cutree(hward, Kfixe)
hcomplete <- hclust(d, method = "complete")
clcomplete <- cutree(hcomplete, Kfixe)
adjustedRandIndex(clward, clcomplete)[1] 0.7540837
clcomplete
clward 1 2 3
1 43 5 0
2 0 47 11
3 0 0 72
fviz_dend(hward, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
show_labels = FALSE, labels_track_height = 0.8)fviz_dend(hcomplete, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
show_labels = FALSE, labels_track_height = 0.8)Question 2
Enoncé
On se concentre maintenant sur la classification hiérarchique par la méthode de Ward. Tracez la courbe des poids de l’arbre et sélectionnez un nombre de classes.
ggplot(data.frame(K = 1:20, height = sort(hward$height, decreasing = T)[1:20]), aes(x = K,
y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") + theme(plot.title = element_text(size = 9))Visualisez la classification obtenue et comparez-la à la répartition des cépages.
Correction
ggplot(data.frame(K = 1:20, height = sort(hward$height, decreasing = T)[1:20]), aes(x = K,
y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") + theme(plot.title = element_text(size = 9)) clward4
1 2 3 4
1 26 20 13 0
2 2 0 18 51
3 0 0 27 21
[1] 0.2818747
fviz_dend(hward, k = 4, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
show_labels = FALSE, labels_track_height = 0.8)Question 3
Enoncé
Mettez en place une classification hiérarchique à partir des coordonnées de l’ACP. Comparez à la répartition par cépages.
Correction
ggplot(data.frame(K = 1:20, height = sort(hwardacp$height, decreasing = T)[1:20]),
aes(x = K, y = height)) + geom_point() + geom_line() + ylab("") + ggtitle("height") +
theme(plot.title = element_text(size = 9))fviz_dend(hwardacp, k = 3, rect = TRUE, rect_fill = TRUE, palette = "npg", rect_border = "npg",
show_labels = FALSE, labels_track_height = 0.8) clwardacp
1 2 3
1 59 0 0
2 6 8 57
3 0 47 1
[1] 0.7605122
Modèles de mélanges
On va maintenant s’intéresser à l’utilisation de modèles de mélanges gaussiens pour déterminer une classification des vins.
Question 1
Enoncé
A l’aide des fonctions mclustBIC() et Mclust(), déterminez une classification des données en considérant
- toutes les formes de mélanges gaussiens
- un nombre de classes maximal à 9 utilisant
- le critère BIC pour la sélection de modèle
Correction
Best BIC values:
VVE,3 EVE,4 VVE,4
BIC -6849.387 -6873.64376 -6885.51547
BIC diff 0.000 -24.25637 -36.12808
On récupère le modèle sélectionné :
On visualise les valeurs du critère BIC pour les différentes formes et les différentes valeurs de \(K\)
Question 2
Enoncé
A l’aide des commandes suivantes, étudiez les probabilités conditionnelles d’appartenance.
Question 3
Enoncé
Visualisez la classification obtenue à l’aide de la fonction fviz_mclust()et comparez-la avec la répartition par cépages.
Question 4
Enoncé
Reprenez les questions précédentes pour la sélection de modèles avec le critère ICL.